clear all; clc;

datafile = 'case300.m'; 
sim = 1;               % 0: reuse older data on initial shedding time
endTime = 1.8;

% Run initial.m & Ymatrix.m
initialize( datafile );
Ymatrix();             % Y3fun: function to change Y by cutting rate

% Simulate before cutting time
simulate(sim, endTime);

%% Get results
global dataSim postSimFunc        % simulation results
global prob nG                      % solution of "case$x.m"
% Set ideal uG and uL; Set No. of generator to display
% sol.uG = []; sol.uL = [];
idx = (1 : length(uG0))'; nidx = length(idx);            % --- choose "delta" to display
% Whole time without cutting gen & load
deltaRec = []; load(dataSim, 'tRec', 'deltaRec', 'coiRec');
cv = deltaRec(:,idx) - repmat(coiRec, 1, nidx); % curve
% Solution of NLP: time phase 3

% Post simulation varification
[tP, dltP, coiP] = feval(postSimFunc, prob.uG0, prob.uL0);
cvP = dltP(:,idx) - repmat(coiP, 1, nidx);
% Plot and Print
plot(tRec, cv, tP, cvP);
axis([0, t3(end), -150, 150]);